library(Mediana)

# Survival-type endpoint (with censoring)
# Survival-type study is an event-driven trial in the presence of censoring.
# A 'Design' object should be added in the data model.

# 1. Data model ----
## The pts are assumed to be recruited at a uniform rate.
## The dropout distribution is exponential.
## Here, dropout means meeting the event (event happened).

## The endpoint type is set to 'event' in the 'OutcomeDist' object:

### Number of events parameters
event_count = c(160, 134)      # Control, Test
randomization.ratio = c(1, 1)  # Control, Test

### Outcome parameters
median_time.control = 14
rate.control = log(2) / median_time.control  # Hazard rate (h)
# h = ln(2) / median_surviving_time 
#   = -ln(surviving_rate) / T0 
#   = -ln(1 - mortality_rate) / T0
outcome.control = parameters(rate = rate.control)
median_time.test = 22
rate.test = log(2) / median_time.test
outcome.test = parameters(rate = rate.test)

### Dropout parameters
dropout.par = parameters(rate = 0.0115)

### Data model
trial.dm = DataModel() +
  OutcomeDist(outcome.dist = "ExpoDist", outcome.type = "event") +
  Event(n.events = event_count, rando.ratio = randomization.ratio) +
  Design(enroll.period = 9, 
         study.duration = 36,
         enroll.dist = "UniformDist",
         dropout.dist = "ExpoDist", dropout.dist.par = dropout.par) +
  Sample(id = "Control", outcome.par = parameters(outcome.control)) +
  Sample(id = "Test", outcome.par = parameters(outcome.test))

# 2. Analysis model ----
## The number of patients and events will be estimated in the study.
## This statistic needs to be specified in a 'Statistic' object.
trial.am = AnalysisModel() +
  Test(id = "Control vs Test",
       samples = samples("Control", "Test"),
       method = "LogrankTest") +
  Statistic(id = "Events Control",
            samples = samples("Control"), method = "EventCountStat") +
  Statistic(id = "Events Test",
            samples = samples("Test"), method = "EventCountStat") +
  Statistic(id = "Patients Control",
            samples = samples("Control"), method = "PatientCountStat") +
  Statistic(id = "Patients Test",
            samples = samples("Test"), method = "PatientCountStat")

# 3. Evaluation model ----
## To compute the average values of the two statistics 
##    (PatientCountStat and EventCountStat), two 'Criterion' are specified.
## And additional Criterion object will be defined to get marginal power.
trial.em = EvaluationModel() +
  Criterion(id = "Marginal Power", method = "MarginalPower",
            tests = tests("Control vs Test"), 
            labels = c("Control vs Test"),
            par = parameters(alpha = 0.025)) +
  Criterion(id = "Mean Events", method = "MeanSumm",
            statistics = statistics("Events Control", "Events Test"),
            labels = c("Mean Events Control", "Mean Events Test")) +
  Criterion(id = "Mean Patients", method = "MeanSumm",
            statistics = statistics("Patients Control", "Patients Test"),
            labels = c("Mean Patients Control", "Mean Patients Test"))

# 4. Clinical Scenario Evaluation ----
## Simulation parameters
trial.sp = SimParameters(n.sims = 5, proc.load = "full", seed = 4293)

## Perform CSE
trial.rst = CSE(trial.dm, trial.am, trial.em, trial.sp)

# 5. Reporting ----
# 5.1 Console
trial.rst.summary <- summary(trial.rst)

# 5.2 Presentation model
# survtrial.pm = PresentationModel() +
#   Project(username = "[Mediana's User]",
#           title = "Survial-type trial",
#           description = "Trial in pts with CRC") +
#   Section(by = "outcome.parameter") +
#   Table(by = "event")

## Generate report
# GenerateReport(presentation.model = survtrial.pm,
#                cse.results = survtrial.rst,
#                report.filename = "SurvivalTypeEndpointTrial.docx")

# 6. Get the simulation data ----
trial.data = DataStack(data.model = trial.dm,
                       sim.parameters = trial.sp)

# i th simulation, j th scenario, k sample
# trial.data$data.set[[i]]$data.scenario[[j]]$sample[[k]]
ExtractSimData <- function(trial.data, scenario, run) {
  # Extract simulation results data (data.set)
  sim_data <- ExtractDataStack(data.stack = trial.data,
                               data.scenario = scenario, 
                               # sample.id = "Control",
                               simulation.run = run
  )$data.set[[1]]$data.scenario[[1]]
  
  data1 <- sim_data$sample[[1]]$data$outcome
  data2 <- sim_data$sample[[2]]$data$outcome
  # exchange 0 and 1 (0 -> 1,, 1 -> 0)
  event1 <- 1 - sim_data$sample[[1]]$data$patient.censor.indicator  
  event2 <- 1 - sim_data$sample[[2]]$data$patient.censor.indicator
  id1 <- rep(sim_data$sample[[1]]$id, length(data1))
  id2 <- rep(sim_data$sample[[2]]$id, length(data2))
  
  rst <- data.frame(value = c(data1, data2), id = c(id1, id2),
                    event = c(event1, event2))
  
  return(rst)
}

sdata <- ExtractSimData(trial.data, 1, 5)
